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We show that transient supersolid quantum states of Rydberg-excitations can be created dynam- 
ically from a Mott insulator of ground state atoms in a 2D optical-lattices by irradiating it with 
short laser pulses. The structure of these supersolids is tunable via the choice of laser parameters. 
We calculate first, second and fourth order correlation functions as well as the pressure to charac- 
terize the supersolid states. Our study is based on the development of a general theoretical tool for 
obtaining the dynamics of strongly interacting quantum systems whose initial state is accurately 
known. We show that this method allows to accurately approximate the evolution of quantum 
systems analytically with a number of operations growing polynomially. 
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The quality of control over atomic systems in state of 
the art experiments is such that one can now address 
and observe quantum evolutions of individual atoms in 
optical-lattices [TJ|2]. Additionally, coherent inter-atomic 
and light-matter interactions can be made strong enough 
to occur on short time-scales compared to incoherent pro- 
cesses. This reveals the system's unitary evolution at the 
individual constituent level which is of fundamental in- 
terest in the study of many-body quantum phenomena. 
Several applications like quantum simulation and quan- 
tum computing schemes also rely on this information 
[3]. The relatively new interest in such non-equilibrium 
dynamics of many-body systems poses serious theoreti- 
cal difficulties due mainly to the exponentially growing 
Hilbert space. This in turn hampers advances in the un- 
derstanding of many-body phenomena such as the elu- 
sive supersolids. This new phase of matter was first sug- 
gested to exist in Helium as the simultaneous existence of 
both diagonal and off-diagonal long-range order [4]. Re- 
cent studies have shown that they might be obtainable 
in optical-lattices [SJ |6] . 

In this letter we propose the transient creation of su- 
persolid quantum states of Rydberg excited atoms from a 
Mott insulator in an optical lattice. These are formed by 
strong laser driving in the presence of long-range dipole- 
dipole interaction. The interaction causes the excitation 
probability of an atom to either be inhibited (blockade) 
or be enhanced (antiblockade) depending on the presence 
of a nearby Rydberg-excited atom [7]. Previously quan- 
tum computation and simulation schemes using similar 
effects have been proposed [5UT2]. Subsequent studies 
have shown that crystal-like dispositions of the Rydberg 
excitations could form the ground states of certain ID 
and 2D lattices [T3HT5] and could also be created dy- 
namically [T6HT9] . Yet these studies have not found su- 
persolidity and have been confined to specific parameter 
regimes. Consequently, the behavior of the system in the 
general case remains largely unknown and thus we here 



develop a novel tool for exploring driven evolutions in 
strongly interacting many-body systems. This is based 
on summing walks performed by an arbitrarily chosen 
subset of the system. In the following we call this method 
walk-sums (WS). 

For a system with N constituents each with d internal 
levels the number of operations involved in computing 
the evolution-operator U scales as d 3N . This results from 
first the number of matrix-elements of U which scales as 
d 2N , and second the number of operations required to ob- 
tain any one of these with a given accuracy, which scales 
as d N . With WS it is possible to approximate analyti- 
cally any chosen pieces of U without any prior knowledge 
about the system under study and with a polynomial 
number of operations in N . The WS method thus solves 
the problem of the second exponential scaling while it 
also bypasses the first one by generating independently 
any desired piece of U. For physical quantities requir- 
ing an exponentially large number of matrix-elements 
of U to be calculated, WS can provide approximations 
by computing randomly chosen pieces of U or only the 
most relevant ones. Therefore we expect WS to work 
well for gapped systems where some configurations are 
unlikely to be populated and can be neglected. WS are 
also more precise than estimates based on truncations 
of the Hilbert-space as it can take into account the ef- 
fect of virtual transitions through configurations outside 
of the truncated Hilbert-space. Additionally, the number 
of operations required per element of U remains exponen- 
tially better with WS than that of a truncated evolution- 
operator. Finally, WS provides a reliable way of getting 
the probability amplitudes of rare events typically inac- 
cessible to Monte-Carlo methods. One can indeed evalu- 
ate the dynamics of a specific piece of the wave function 
and concentrate the computational effort on obtaining a 
very high accuracy on this single piece. The reasoning we 
present here is a special case of a very general procedure 
for working out elements of general matrix functions by 



2 



summing paths on a graph [20] . 

WS are based on splitting a many-body system into a 
set of constituents S' whose dynamics is frozen and then 
computing the evolution of the remaining few particles 
S. The surrounding S' being perfectly known, all inter- 
actions with S can be exactly evaluated and S evolves 
through a small Hamiltonian depending on the configu- 
ration of S' . By expressing the true many-body dynamics 
in terms of such simple situations with a frozen S' and a 
few evolving constituents S one can solve the Schrodinger 
equation for large systems. This might seem similar in 
essence to mean field theory where an atom of interest 
interacts with a field resulting from the mean behavior of 
all other particles. The difference is that here we make 
the mapping from many-body to few-body dynamics ex- 
act, that is we make S interact with all possible fields it 
could be subjected to depending on the configuration of 
S' . This will effectively perform some average and mean 
field behaviors will be recovered but fluctuations around 
this mean will also be present. 

We project the system onto a specific configuration 
of all constituents in S f by applying the operator £ M = 
®2s called a projector-lattice. Here /i denotes a 
basis state configuration in S f and the identity is applied 
to S. This operator satisfies the closure relation ^2 = 
X which we insert into the system evolution operator J7, 
written as a product of infinitesimally small steps in time 
St. This leads to U(t) = E M S V lim ^o ]ir{E M ^1^' 
with SU = U(St) and time t = mSt. Expanding this 
product yields terms like 

iyUit) = lim e v SUi v SU . . . e v SU + 

lim e u SU . . . i u SUe u SU . . . eJU + . . . (1) 

The first term in this expansion evolves the system from 
to St, then e v projects S 1 onto |z/), followed by evo- 
lution for St, etc. This freezes S f by continuous (Zeno) 
measurement of e v in the limit St —> 0, while S evolves 
freely. The other terms in this expansion describe any 
number of consecutive Zeno measurements of different 
configurations of S' switching at all possible times. For 
instance, the second term of Eq. corresponds to one 
change from configuration \i to v. Thus, in this expan- 
sion we consider S to evolve in an environment S' which 
evolves stroboscopically and simultaneously through all 
possible configurations. 

Provided that S f evolves from an initial configuration 
\i to a final configuration v we thus write the evolution 
operator for sub-system S as e v U(t)£^ = U v ^^{t)^\v){ii\ 
with (ft 1) 

W n (G) J ° J ° J ° 

H^ Vn _ ie - iH ^-i . . . H m ^e- iH ^dh . . . dt n . (2) 



We note that this expression only contains d x d ma- 
trices. The index n indicates the number of jumps 
undergone by S' between and t and the sum over 
W n (G) = {rjo = /i, 7/1 • • • Vn-iiVn = ^} contains all pos- 
sible strings of n consecutive jumps starting at configu- 
ration ji and ending at v. The integrals represent con- 
tinuous sums over all the possible jumping times with 
one integral per jump. The d x d matrices H Vj are ef- 
fective Hamiltonians evolving S for a given configura- 
tion r]j of S' while the matrices H f]j ^ r]j _ 1 describe the 
effect of a jump on S. They are given by e^.Re^. = 
\Vj)(Vj\ irj.Hi^ = \ri j )(ri j - 1 \®H rij ^ ri ._ 1 . The 

operator U v ^^(t) is thus the conditional-evolution op- 
erator for S knowing that initially S' was in state 
and in state \v) at time t. The time integrals of Eq.([2| 
are convolutions and as a consequence the expression 
of conditional-evolution operators in the Fourier domain 
only involves additions and multiplications of d x d ma- 
trices M Vj {uj) = JT[#(t)e~ ijH V], with 0(t) the Heaviside 
function and TT the Fourier transform. In the Fourier 
domain Eq.Q becomes 

[/■„«_„ = $3 r n M v H v ^ Vn _ 1 M Vn _ l ...H v ^ ll M ll . 

« W n (G) 

(3) 

In the case of 2 x 2 matrices, we find the M Vj to have a 
simple universal expression in terms of the H r]j and that 
all matrix elements of any product of these M Vj are ra- 
tios of polynomials with analytically known roots. It fol- 
lows that we can always analytically perform the inverse 
Fourier transform back into the time domain for 2 x 2 ma- 
trices. We will see that this formulation of the dynamics 
is efficient, but a trade-off is that we obtain exponentially 
little information. Indeed, we really compute s u U(t)s^ 
i.e. only d 2 out of d 2N elements of the full evolution- 
operator U. Computing the conditional-evolution opera- 
tors for a large number of final configurations is possible 
because we have analytic expressions, and we can there- 
fore approximate a large number of pieces of U(t). We 
will follow this approach for Rydberg excitations below. 

While every element of Eq. (|3| can be worked out ef- 
ficiently it will in general be difficult to distinguish all 
the possible strings of jumps. We thus proceed by map- 
ping the possible strings of jumps W n (G) to walks on 
a graph G. We construct a graph G as follows (i) for 
each configuration rjj we draw a vertex v Vj , and (ii) for 
each H r]j ^ r]i / we draw an edge between vertices v Vj 
and v Vi . Now all the possible successions of jumps be- 
tween the initial and final configurations of S' correspond 
to all the possible walks W n (G) on G between the ini- 
tial and final vertices and v v . The length of a walk 
is equivalent to the number of jumps in S' . We note 
that it is possible to derive Eq. Q directly from the 
power-series expansion of the matrix-exponential of the 
Hamiltonian by exactly summing terms of the form H™. 
which appear as loops on the vertices of G [25] . The re- 
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FIG. 1: Walks of order 1 and 3 contributing to £/j<_o(£), 
the conditional-evolution operator between the initial con- 
figuration with no excitations £ = (vertex 0) and the fi- 
nal configuration with one excitation £ = 1 on atom j (ver- 
tex rj). The corresponding mathematical operations are 
£Vo MjH^oMo + ^{M^H^oMoHo^iMkH^oMo + 
MjH^Mj^^iMjH^oMo} + . . . 

moval of these loops makes our expansion different from 
a pure power-series expansion. It leads to a significant 
speed-up when calculating individual elements U v ^^ and 
guarantees that a truncation of Eq. Q at order K is at 
least as accurate as a similar truncation of the power- 
series. To make this statement quantitative, we com- 
pare the number of floating point operations required 
per matrix element g with that of the power-series ex- 
pression of U . We obtain the number of walks of a given 
length K between two vertices of a graph from its ad- 
jacency matrix Aq [21] . The contribution to Eq. pi of 
order K requires the multiplication of (2K + 1)(v\Aq\ih) 
matrices with (v\Aq\/jl) the number of walks of length 

K which is upper bounded by Yll (jz^) — ^ ' ^ or 
most physically relevant Hamiltonians we found q = 1,2. 
Therefore g < d~ 2 (2Kd 3 + d 2 ) {q (d-l) q } K , where 

(2Kd 3 + d 2 ) is the number of operations required to mul- 
tiply the 2K + 1 matrices of a walk of length K and to 
add the result to the other walks. We compare this to 
the corresponding value of the Taylor-series expansion 
gT = (K — l)d 3N / d 2N c± Kd N , and since g is polynomial 
in A/", the ratio g/gr —> with increasing size N — >• oo. 



fast laser pulses to several /is but can safely be neglected 
on the much shorter time scales considered here [22]. In 
this limit the atoms are described by the Hamiltonian 

H = J2i^Pi - § (Ti + T}) + AijPiPj}, (5) 

i j^i 

with A the laser detuning, Q the Rabi-frequency, Ti = 
\g)i(r\ and Pi = \r)i(r\. We choose S to be the atom at 
the center of the lattice and construct a graph G where 
each vertex represents a configuration with £ Rydberg 
excitations. The Hamiltonian drives transitions £ —> £±l 
and so G is a linear graph with the vertex representing no 
excitations at its end. All Hi±i^£ = — (Q/2)T S while the 
Hamiltonian with £ excitations at positions j = {ji . . . ji} 
in S' is given by 

^=^ + E^) + (_n/ 2 a + eXj)- (6) 

Figure [l] shows an example of how we sum walks on G to 
determine elements of U. We calculate all conditional- 
evolution operators whose final configuration has up to 
six simultaneous excitations anywhere in the lattice. This 
gives analytical approximations to 2 2 x Y^ z =o (^) e ^ e ~ 
ments of the many-body evolution operator U (t) result- 
ing from billions of walks on the graph. For N > 100 we 
limit our calculations to three excitations and randomly 
chosen subsets of configurations with four or more excita- 
tions. We typically obtain the many-body wave function 
\ij){t)) for 2D lattices with A^ « 6600 with moderate com- 
putational effort [26 for arbitrary sets of parameters in 
the Hamiltonian. We then calculate the order-parameter 
(Tj) and various correlation functions and check their 
convergence by varying the randomly chosen samples of 
configurations with 4 < £ < 6. Contributions with i > 6 
are not considered. 



We now apply WS to lattices of strongly interacting 
Rydberg atoms. Initially the atoms are assumed to be in 
a pure Mott insulating state of the form \gg . . . g). The 
atoms are arranged in a regular 2D pattern by trapping 
them in a deep optical lattice with one atom in inter- 
nal ground state \g) occupying each lattice site. A laser 
drives the atoms to highly excited Rydberg states \r) 
which are strongly interacting over long distances via a 
dipole-dipole interaction of the form [3 [8] 

A i:j = ^Treoi^O" 1 ^^ ( 4 ) 

with \±i the dipole moment of atom i and Rij the relative 
distance between atoms i and j. The laser and dipole- 
dipole interactions we consider are in the MHz-GHz range 
and induce dynamics fast compared to incoherent pro- 
cesses and the motion of the atoms. These will limit the 
lifetime of Rydberg excited quantum states created by 



A simple picture of the density-density correlation 
function ^(s, j) = (r - 3 r 'j) / '(r s ) (r j) is available at the low- 
est order, where we find it to be given by 



n 2 



- sin 



2 n 2 

x 2 



- sin 



{ X 2 sm[^f ty 



(7) 

with x 2 = ^ 2 + A 2 and x % = ^ 2 + (A + A sj ) 2 . This can 
be understood from conditional-probabilities ^(s, j) = 
(1 - (r s /rj) + {r s /gj)){rs/rj)/{r s /gj) with e.g. (r s / gj ) 
being the probability of finding atom S excited provided 
j is in the ground state. At this low order no atom is ex- 
cited in the lattice except for S and j and the conditional- 
probabilities are expected to be identical to those of the 
two- atom problem, just as found here. Thus as long as 
the probability of having more than one excitation in the 
vicinity of S can be neglected, higher orders should be ir- 
relevant. This result is accurate for up to ~ 27r-pulses in 
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FIG. 2: (Color online) Correlation-functions multiplied by the number of pairs a) g[(s, j) = (T^T? +T s Tj) - (T s + T s f )(Tj + T/), 
b) pi(s, j) = (TsTf + T 8 + Tj) - (T s + T^T- + Tj 1 ") and c) p 2 (s, j) over a JV « 6600 atoms square lattice with lattice spacing 
L = 1.5/im, one atom/pixel, s is the central one. #2(s, j) ~ except at the few sites where s and j form a free pair, d) Locations 
of the dominant peaks of g± — {r s T^r\^r\) / (r s ) (rj) (rk) (n) (black squares) and #4 = (g s T^r\^r\) / (g s ) (rj) (rk) (n) (red disks) for all 
j, k forming free pairs with s. Parameters : Rydberg-state principal number n — 40, £1 = 30MHz, A = 0, Qt = 8tv, ~ 0.437T 
and = 7r/2, convergence indicates a precision ~ 10 -4 . 



strongly blockaded situations and provides similar values 
as the ^-expansion for ID systems [23]. Because Eq.Q 
is just a two-body picture however, it completely fails to 
account for correlations between atoms that do not inter- 
act directly and where thus four body processes become 
dominant. 

This situation is e.g. realized in square lattices when 
the laser angles with respect to the lattice plane (0, (j)) 
fulfil 

3sin 2 (0)(;fccos(0) + j^sin^)) 2 = j 2 + j 2 , (8) 

where j = (j x ,jy) are the lattice coordinates of atom 
j and s = (0,0). Then the interaction term A s j = 
and in general \A s y\ < ft on the line joining s and j. 
From now on pairs s, j fulfilling Eq.([8| will be called free 
pairs. In the presence of such pairs, we observe states 
with non-zero order-parameter (Tj). As shown in Fig. 2 
they also display substantial and correlated ^(sj) and 
#i(s,j), anticorrelated with <7i(s,j) and all peaking pre- 
cisely where s and j form a free pair. Additionally, we 
find (r s ) = Xlj( r s r j) within our numerical accuracy in- 
dicating that excitations occur at least in pairs, and this 
sum being dominated by a 99% contribution from free 
pairs. By tuning and (j) one can choose which pairs of 
atoms are free and hence the structure of the correlations. 

To further determine the nature of the observed states 
we compute their pressure p^D = — (dE/dA)T,Nf R where 
E and A = NL 2 are the energy and area of the system, 
respectively. The derivative is taken at constant tem- 
perature T = and number of excitations N/r, with 
fn = EMto)\P^(t ))/N the Rydberg fraction. To 
obtain P2D we turn off the laser at time to and con- 
sider a contraction of the system L(t). The resulting 
Hamiltonian H 1 commutes with itself at any time as 
well as with any of the P%Pj ^ nus leaving fn unchanged 
while E(t) = Eij^uW(^(*o)|^j|^(*o)> changes only 
through L(t). This yields P2D = (3/2)E/A. As discussed 
above the sum in the energy is dominated by free pairs 



for which A^ = and we thus find these states to have 
a very small pressure. Physically this is because Eq.(|8| 
is independent of the lattice spacing L. Consequently, 
during a contraction, atoms forming free pairs never in- 
teract directly. We expect this small pressure to increase 

— 1/2 

the lifetime of the observed states r ex p 2B compared 
to states containing interacting excitations. 

The presence of both diagonal and off-diagonal long- 
range order together with a nearly vanishing quantum 
pressure qualifies the observed states as supersolids [24] . 
Their translational symmetry can only become appar- 
ent from n > 4-body correlation-functions as shown in 
Fig.j2]-d). Finally we note that around = and 7r/2, 
Eq.([8f is a robust condition as errors of Acj) and A# 
on the laser angles produce a small residual interaction 
I Ay I < 2V^/ii/ij(47re i?? j )" 1 A# in free pairs. 

In this article we have analysed many-body Rydberg- 
excited quantum states created by short laser pulses from 
an atomic Mott insulator of ultracold atoms. We have 
calculated their pressure, and first, second and fourth or- 
der correlation functions to show that these states ful- 
fil the properties of a supersolid for specific laser se- 
tups. Our calculations were based on WS which provide 
a versatile tool for investigating the coherent dynamics of 
many body quantum systems whose initial wave function 
is known. WS are part of a generic method to evalu- 
ate general matrix functions and can e.g. be extended 
to analyse continuous-time isotropic quantum random 
walks on any lattice geometry and be shown to reduce to 
path-integrals for systems with continuous-only degrees 
of freedom. 
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